clear
version 14
set more off
cap log close

mata: mata clear
mat drop _all

gl OUT 	"${LOCDRIVE}\Results"
gl TEMP	"${LOCDRIVE}\Tempdata"
gl FIG "${LOCDRIVE}\Figures"

loc r1_hh hhid
loc r1_weight hh_weight
loc r1_indiv sbmemno
loc r2_hh y2_hhid
loc r2_weight y2_weight
loc r2_indiv indidy2
loc r3_hh y3_hhid
loc r3_weight y3_weight
loc r3_indiv indidy3

loc ea ea
loc strata strataid
loc admin0 macroregion
loc admin1 region 		
loc admin2 district  	
loc admin3 ward

loc numeraire other

loc cre 1 /*include correlated random effects*/
loc crelev hh_c /*level of correlated random effects*/

loc servicelab "Services"
loc foodlab "Food"
loc goodlab "Goods"
loc otherlab "Other"

loc 3g1 food
loc 3g2 good
loc 3g3 service 

	
log using "${TEMP}\Analaysis_compfigs", replace


/*#Set which specifications to compare for the figures*/

loc  mainspec 	SL_food1good1service1_np5 
loc nf_mainspec = 1 
loc mainspec_lab "IV Estimates"
loc  altspec 	SL_food0good0service0_np5 
loc nf_altspec = 1
loc altspec_lab "IVs as controls"

loc S_mainspec = `nf_mainspec' + 2
loc S_altspec = `nf_altspec' + 2
loc S_max = max(`S_mainspec',`S_altspec')

loc nbins = 10


** Set up paramlist

loc dsex_paramlist "dsf_ex"
loc dsex_collapselist_sd "dsf_ex_se=dsf_ex"
loc eeex_paramlist "eef_ex"
loc eeex_collapselist_sd "eef_ex_se=eef_ex"
loc eeepr_paramlist ""
loc eeepr_collapselist_sd ""

forvalues s=1/`S_max' {   
	loc dsex_paramlist `dsex_paramlist' ds`s'_ex
	loc eeex_paramlist `eeex_paramlist' ee`s'_ex
	
	forvalues k=1/`S_max' {
		loc eeepr_paramlist `eeepr_paramlist' ee`s'_epr`k'
		}		
	}

display "`dsex_paramlist'"
display "`eeex_paramlist'"
display "`eeepr_paramlist'"

loc paramlist `dsex_paramlist' `eeex_paramlist' `eeepr_paramlist' 

	
** Elasticity figures by quintile

use "${OUT}\reg3_elasticities_`mainspec'", clear
keep if round==3 /*This is new*/
g spec = "`mainspec'"
tempfile mainspec_r3only
save `mainspec_r3only', replace

use "${OUT}\reg3_elasticities_`altspec'", clear
keep if round==3 /*This is new*/
g spec = "`altspec'"
tempfile altspec_r3only
save `altspec_r3only', replace

use `mainspec_r3only', clear
append using `altspec_r3only'

tempfile reg3_elasticities_bin
save `reg3_elasticities_bin'

foreach pop in urb rur all {
	foreach param in `paramlist' {
		use `reg3_elasticities_bin', clear

		keep if `pop'==1
		xtile bin = y1, n(`nbins')

		sort bin spec `param'
		
		cap su `param'
		if _rc==111 {
			g `param' = .
			g `param'_se = .
 			}
		
		collapse (median) `param' `param'_se, by(bin spec)

		g pop = "`pop'"
		keep bin pop spec `param' `param'_se

		tempfile `param'_medval_`pop'
		save ``param'_medval_`pop'', replace		
		}
	}
	
	
clear
set obs `nbins'
g bin = _n
expand 3

g count = 1
bysort bin : g popc = sum(count) 
g pop = "all" if popc==1
replace pop = "urb" if popc==2
replace pop = "rur" if popc==3
replace count = 1

expand 2
bysort bin popc : g specc = sum(count)
g spec = "`mainspec'" if specc==1
replace spec = "`altspec'" if specc==2

foreach S in 3 6 {
	loc shift = 1
	forvalues s = 1/`S' {
		g order`S'_`s'm = (bin-1)*(`S'*2+1) + `shift' if spec=="`mainspec'"
		loc shift = `shift' + 1
		g order`S'_`s'a = (bin-1)*(`S'*2+1) + `shift' if spec=="`altspec'"
		loc shift = `shift' + 1
		}
	}
	
drop popc count

foreach pop in urb rur all { 
	foreach param in `paramlist' {
		foreach g in 1 2 3 {
			merge 1:1 bin pop spec using ``param'_medval_`pop'', nogen update
			}
		}
	}
	
foreach param in `paramlist' {
	g `param'_sim_lo = `param' - `param'_se*1.96
	g `param'_sim_hi = `param' + `param'_se*1.96		
	}

save "${OUT}\elasticity_summary_byq_comp_`mainspec'_`altspec'", replace


** Set up the Labels margins and orders 

loc dslab "Budget share semi-elasticity"
loc dsyline 0
loc ymin_ds -0.3
loc ymax_ds 0.35
loc yskip_ds 0.1
loc eelab "Expenditure elasticity"
loc eeyline 1
loc ymin_ee 0
loc ymax_ee 3
loc yskip_ee 0.5
loc eeeprlab "Hicksian price elasticity"
loc eeepryline 0
loc ymin_ee_epr -3
loc ymax_ee_epr 3
loc yskip_ee_epr 0.5


foreach S in 3 6 {
	loc xlab`S' ""
	loc xlinecode`S' ""
	loc o = ((`S'*2)+1)/2
	forvalues b=1/`nbins' {
		loc xlineloc = `o' + ((`S'*2)+1)/2
		loc xlinecode`S' `xlinecode`S'' xline(`xlineloc', lwidth(vthin) lpattern(dash) lcolor(gray))
		loc xlab`S' `xlab`S'' `o' "bin `b'"
		loc o = `o'+(`S'*2)+1
		}
	display `"The label for S=`S' is: `xlab`S''"'
	}

foreach param in ds ee {
	foreach j in s g {
		g `param'`j'_ex = .
		g `param'`j'_ex_sim_hi = .
		g `param'`j'_ex_sim_lo = .
		}
	}
		
foreach spec in mainspec altspec {
	foreach param in ds ee {
		replace `param'f_ex = `param'1_ex if spec=="``spec''"
		replace `param'f_ex_sim_hi = `param'1_ex_sim_hi if spec=="``spec''"
		replace `param'f_ex_sim_lo = `param'1_ex_sim_lo if spec=="``spec''"
		replace `param'g_ex = `param'2_ex if spec=="``spec''"
		replace `param'g_ex_sim_hi = `param'2_ex_sim_hi if spec=="``spec''"
		replace `param'g_ex_sim_lo = `param'2_ex_sim_lo if spec=="``spec''"
		replace `param's_ex = `param'3_ex if spec=="``spec''"
		replace `param's_ex_sim_hi = `param'3_ex_sim_hi if spec=="``spec''"
		replace `param's_ex_sim_lo = `param'3_ex_sim_lo if spec=="``spec''"
		}		
	}
		
	
** Expenditure Elasticities 

foreach param in ds ee {
	foreach pop in urb rur all {
		graph twoway 	(rcap `param'f_ex_sim_hi `param'f_ex_sim_lo order3_1m, color(red%15)) ///
							(rcap `param'f_ex_sim_hi `param'f_ex_sim_lo order3_1a, color(red%15)) ///
							(scatter `param'f_ex order3_1m, color(red%15) msymbol(O))  ///
							(scatter `param'f_ex order3_1a, color(red%15) msymbol(D))  ///
							(rcap `param'g_ex_sim_hi `param'g_ex_sim_lo order3_2m, color(green%50)) ///
							(rcap `param'g_ex_sim_hi `param'g_ex_sim_lo order3_2a, color(green%50)) ///
							(scatter `param'g_ex order3_2m, color(green%50) msymbol(O)) ///
							(scatter `param'g_ex order3_2a, color(green%50) msymbol(D)) ///
							(rcap `param's_ex_sim_hi `param's_ex_sim_lo order3_3m, color(blue%80)) ///
							(rcap `param's_ex_sim_hi `param's_ex_sim_lo order3_3a, color(blue%80)) ///		
							(scatter `param's_ex order3_3m, color(blue%80) msymbol(O)) ///
							(scatter `param's_ex order3_3a, color(blue%80) msymbol(D)) ///
							if pop=="`pop'" ///
							, legend(order(- "`mainspec_lab':" 3 "Food" 7 "Goods" 11 "Services" - "`altspec_lab':" 4 "Food" 8 "Goods" 12 "Services") rows(2)) title("``param'lab' w.r.t. total exp., `pop' pop") xlabel(`xlab3', noticks) xtitle("Log real expenditures") yline(``param'yline') ylabel(`ymin_`param'' (`yskip_`param'') `ymax_`param'', nogrid) yscale(range(`ymin_`param'' `ymax_`param'')) `xlinecode3'
			graph export "${FIG}\Elast_byq_foodagg_`param'ee_`pop'_comp_`mainspec'_`altspec'.pdf", replace 
			}
		}
			
			
** Price Elasticities 
	
foreach param in ee {
	foreach pop in urb rur all {
		forvalues g = 1/`S_max' {
			graph twoway 	(rcap `param'1_epr`g'_sim_hi `param'1_epr`g'_sim_lo order3_1m, color(red%10)) ///
							(rcap `param'1_epr`g'_sim_hi `param'1_epr`g'_sim_lo order3_1a, color(red%10)) ///			
							(scatter `param'1_epr`g' order3_1m, color(red%10) msymbol(O))  ///
							(scatter `param'1_epr`g' order3_1a, color(red%10) msymbol(D))  ///
							(rcap `param'2_epr`g'_sim_hi `param'2_epr`g'_sim_lo order3_2m, color(green%50)) ///
							(rcap `param'2_epr`g'_sim_hi `param'2_epr`g'_sim_lo order3_2a, color(green%50)) ///		
							(scatter `param'2_epr`g' order3_2m, color(green%50) msymbol(O)) ///
							(scatter `param'2_epr`g' order3_2a, color(green%50) msymbol(D)) ///
							(rcap `param'3_epr`g'_sim_hi `param'3_epr`g'_sim_lo order3_3m, color(blue%80)) ///
							(rcap `param'3_epr`g'_sim_hi `param'3_epr`g'_sim_lo order3_3a, color(blue%80)) ///
							(scatter `param'3_epr`g' order3_3m, color(blue%80) msymbol(O)) ///
							(scatter `param'3_epr`g' order3_3a, color(blue%80) msymbol(D)) ///
							if pop=="`pop'" ///
							, legend(order(- "`mainspec_lab':" 3 "Food" 7 "Goods" 11 "Services" - "`altspec_lab':" 4 "Food" 8 "Goods" 12 "Services") rows(2)) title("``param'eprlab' w.r.t. `3g`g'' price, `pop' pop") xlabel(`xlab3', noticks) xtitle("Log real expenditures") yline(``param'epryline') ylabel(`ymin_`param'_epr' (`yskip_`param'_epr') `ymax_`param'_epr', nogrid) yscale(range(`ymin_`param'_epr' `ymax_`param'_epr')) `xlinecode3'
			graph export "${FIG}\PriceElast_byq_`param'_epr_`3g`g''_`pop'_comp_`mainspec'_`altspec'.pdf", replace  
			}
		}
	}

	